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1. Introduction 

QCD at high temperature and density has rich phase structure, and the nature of the phase 
transition changes depending on the quark mass and the number of flavors. The quark mass de¬ 
pendence of the QCD phase transition is important not only for the theoretical interest in the phase 
structure of QCD with various quark masses but also for the investigation of finite density QCD or 
QCD near the chiral limit, in which the numerical studies are difficult. From the numerical study 
of (2+l)-flavor QCD, the QCD phase transition is considered to be crossover at low density but 
is expected to turn into a first order transition at high density. Finding the endpoint of first order 
transition is then crucial for establishing the above expectation. Unfortunately, it is still extremely 
difficult to simulate quarks with physical mass at high density. But accessing the endpoint becomes 
easy if one extends QCD in appropriate directions, and importantly if such an extension is smooth, 
one can extract information on original QCD by suitable extrapolations [|I[ |2j]. 

Moreover, the nature of the phase transition in the chiral limit of 2-flavor QCD is a long¬ 
standing problem. The standard expectation is of second order, but the order is not conclusive 
due to the difficulty of the numerical study in the chiral limit. One of the good approaches is to 
investigate the boundary of the first order transition region. If the critical value of the strange quark 
mass does not go to infinity in the up and down quarks massless limit, the transition of the massless 
2-flavor QCD is not of first order. However, it is difficult to study it, because the first order region 
is very small in (2+l)-flavor QCD and simulations with very small mass are required. 

In Ref. [Q], the boundary of the first order region is studied in many-flavor QCD, motivated 
by a feasibility study of the electroweak baryogenesis in technicolor theories constructed by many- 
flavor gauge theory. They studied QCD with two degenerate light quarks of the mass m\ and the 
chemical potential p\ and Af massive quarks with and p^, and found by measuring probability 
distribution function that the critical massive quark mass becomes larger as Af increases. Therefore, 
the investigation of the critical line becomes easier as Nf increases. This extension of QCD is also 
useful for the study of the phase transition of massless 2-flavor QCD. If the critical mass for the 
massive Nf flavors remains finite in the massless limit of two light flavors, it gives a strong support 
for the second order transition. A similar approach has been tried in Ref. [[Jj]. 

This paper consists of two parts. The first part deals with the light quark mass dependence of 
the critical massive quark mass separating the first order and crossover regions in (2 + Nf )-flavor 
QCD. We perform simulations with 2 flavors of improved Wilson quarks. The effect from the 
dynamical Af-flavors are added by a reweighting method assuming that the Af-flavors are heavy. 
We then investigate the critical heavy quark mass as a function of the light quark mass through 
the shape of distribution functions and discuss the nature of phase transitions in the chiral limit of 
2-flavor QCD. 

The second part is the chemical potential dependence. In particular', we extend the real chemi¬ 
cal potential (ju) to complex value. The phase transition of 2-flavor QCD with finite mass at p = 0 
is crossover. But, when the chemical potential is introduced, the shape of the distribution changes 
and the distribution function will become first-order-transition-like. For the complex p, this prop¬ 
erty make a singularity, which is called “Lee-Yang zero” [||]. Lee and Yang proposed the method 
to investigate the nature of phase transitions from the singularities in the complex parameter plane, 
and applications to finite density QCD have been discussed in Refs. M- 
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In the next section, we explain our method to identify the nature of phase transitions via 
the distribution function. We then argue the light quark mass dependence of the endpoint of the 
first order transition in Sec. ||. The singularities in the complex plane are discussed in Sec. Jlj. 
Conclusions are given in Sec. ||. 


2. Critical point by a histogram method 


We study QCD with two degenerate light quarks of the mass m\, the chemical potential f.i\ and 
Nf heavy quarks. We define the probability distribution function of average plaquette value, 


) = JSdUS’y/Sty 8(P — P) 


— f) e Sg 


N { 


{®U 8{P-P) e 6liN ^ p (2.1) 

■' /=i 


where S g and S q are the actions of gauge and quark fields, respectively, and M is the quark matrix. 
For simplification, we adopt M which does not depend on /j explicitly. N s - ac = N 2 x N, is the number 
of sites and /3 = 6/gg is the simulation parameter. P defined by P= —S g / (6/V s i te /3) is 1 x 1 Wilson 
loop for the standard plaquette gauge action. 8(P — P) is the delta function, which constrains the 
operator P to be the value of P. We moreover define the effective potential, V^(P:fi.mf.pf ) = 
-iaw(P\P,m f ,Hf). 

Denoting the potential of 2-flavor at /u= 0 by Vq(P\ j8), that of (2 + W)-flavor is written as 


VeffC P-,P,m f ,n f ) = V 0 (P\Po)-]nR(P\P,m f ,Hf,Po), (2.2) 


with 


ln/?(P;j8,m/,jU/;j8o) 


603 -/3 o )fVsiteP + ln 


/detM(mi,ju)\ 2 -py det M(trif,Hf) \ 

\ detM(mi, 0 ) / detM(oo,Q) / ^ 


where (■ • ^(prnxed) = (8(P — P) ■ ■ ')fJ(8(P — P))p n and (■■■)p 0 means the ensemble average over 
2-flavor configurations generated at /3o, m\ and vanishing jl\. f5 o is the simulation point, which may 
differ from fi in this method. 

At a first order transition point, V e s shows a double-well shape as a function of P, and equiva¬ 
lently the slope of the potential dV e ff/dP shows an S-shape. Since /j appears only in the linear term 
of P on the right hand side of Eq. ( |2.3[ ), the shape of the slope dV e ff/dP is independent of /3, i.e. a 
change of /J just shifts the overall constant of the slope of Feff [§]■ Although /3 must be adjusted to 
the first order transition point to observe the double-well potential, the fine tuning is not necessary 
if we investigate the slope. 

The derivative dVo/dP can be measured easily from the peak position of the plaquette his¬ 
togram [§]. When one performs a simulation at /3o, the slope is vanishing at the minimum of 
Vo(P; /3o), and the value of P at the minimum can be estimated by (P)p lt approximately. Hence, we 
obtain dVo/dP at (P)p ll by 

^((P) ft , j S) = -6( J 8-j3o)tV si te. (2.4) 

We therefore focus on the slope of the effective potential to identify the nature of transitions. 
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Figure 1: Left: Plaquette histogram of 2-flavor QCD wo(P;/3, fe.O) at k = 0.145. Right: In RiP'JiS)) as a 
function of P for h = 0.05 - 0.50. 


3. Light quark mass dependence of the critical heavy quark mass 


In this section, we discuss the light quark mass dependence of the critical hopping parameter of 
heavy quarks terminating the first order transition. It is important to investigate whether the critical 
hopping parameter vanishes at a finite light quark mass or not. If the critical line crosses the line 
of Kh = 0, the 2-flavor chiral limit is in the first order region, where fq, is the hopping parameter 
of heavy-flavors being proportional to the inverse mass. Restricting the calculation to the heavy 
quark region near jq, = 0, the second determinant for Nf flavors in Eq. ( |2.3| ) is approximated by the 
leading order of the hopping parameter expansion, 


In 


detM(ic h , i Uh) 

detAf(0,0) 


= 288Al site J^P + \2N]{2K h ) N < (cosh( J u h /r)^ R + /sinh( J u h /r)f2 I ) + • • • (3.1) 


for the standard Wilson quark action [|1 Q|]. Q R and Lfi are the real and imaginary part of the 
Polyakov loop, respectively. For improved gauge actions such as S g = —6A2 s i te /3 [co (plaquette) + 
ci (rectangle)], additional ci x 0{ fc 4 ) terms must be considered, where ci is the improvement coef¬ 
ficient and co = 1 — 8ci. However, since the improvement term does not affect the physics, we will 
cancel these terms by a shift of the coefficient ci. It is shown in Ref. [j|] that the hopping parameter 
expansion is applicable for large Nf. 

Denoting h = 2Nf(2K\f) N ', we obtain lnR(P;/3, ffh,0;/3o) = \nR(P:h.()) + (plaquette term) 
+0(Kjf' +2 ) for H\ = /r h = 0 with 


R(P-,h,0) = (exp[6/W s 3 D R ]) (p:fixed ft) , (3.2) 

where R(P',h, 0) is given by the Polyakov loop and is independent of /3 q. The plaquette term 
does not contribute to d 2 V e ff/dP 2 and can be absorbed by shifting jS —>• jS* = jS +48A r ffC^ for 
Wilson quarks. Moreover, one can deal with the case with non-degenerate masses by adopting 
h = 2£jli (2Kf) N ‘ for the Wilson quark action or h = (1 /4) Y./L i (2mf)^ N ' for the staggered quark 
action. Thus, the choice of the quark action is not important. In the following, we discuss the mass 
dependence of R through the parameter h. 
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16 3 x4, K=0.145, c m =1.650 16 3 x4. K=0.1475, c m =1.677 



16 3 x4. K=0.150, Csw =1.707 16 3 x4. K=0.1505, c sw =1.712 



Figure 2: dV e ff/dP(P',P,K,h) as a function of P at jq = 0.1450(top left), 0.1475 (top right), 0.1500 (bottom 
left) and 0.1505 (bottom right). 


We perform simulations of 2-flavor QCD with the clover-improved Wilson quark and Iwasaki 
gauge actions 1 . We take four different values of light quark masses ranging from fq = 0.145 to 
0.1505. The corresponding ratio of pseudo-scalar and vector meson masses is mps/my ~ 0.6647, 
0.5761, 0.4677 and 0.4575 for ?q = 0.145, 0.1475, 0.150 and 0.1505, respectively. The data are 
taken at about 30 values of (5 around the pseudo-critical point at each K /, and at each simulation 
point 10,000 to 40,000 trajectories have been accumulated. 

In the calculation of R(P\h,0), we use the delta function approximated by 8(x) ~ I /(A^/ji) 
exp[— (jc/A) 2 ], where A is selected consulting the resolution and the statistical error. Because 
R(P;h, 0) is independent of /3, we mix all data obtained at different /3 as was done in Ref. [j81. 
The histograms of plaquette value are plotted in the left panel of Fig. p]. The results for In R(P; h, 0) 
at K\ = 0.145 are shown in the right panel of Fig. p] for h = 0.05 - 0.50 at interval of 0.05. A rapid 
increase is observed around P ~ 1.6, and the gradient becomes larger as h increases. 

The derivative dlnR/dP is calculated by fitting In R to a / 2 lh -order polynomial of P in an ap¬ 
propriate fit range. Using the equation, 


^Veff 

clP 


dVo 

dP 


d(lnR) 

dP 


+ (const.), 


1 The improvement parameters are almost the same as Ref. [|TT|j 


(3.3) 
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2 

( m i/ m p) 


Figure 3: The critical h as a function of ( m K /m p ) 2 (left) and of wpcac (right). The red liens are lower and 
upper bounds of the critical value of h obtained by dynamical simulations of degenerate 16-flavor QCD. 


we compute dV e ff/dP for each jq. The results for various h are shown in Fig. Q, where dVq/dP 
is computed by Eq. ([2d]). The shape of dV^/dP is independent of because d 2 V c a/dP 2 is /J- 
independent. The first derivative dV e n/dP is the monotonically increasing function of P when h is 
small, indicating that the transition is crossover. However, the shape of dV e s/dP turns into an S- 
shaped function at h ~ 0.3, which means that the system undergoes first order transition. The same 
analysis has been done in Ref. [|j] using the p4-improved staggered fermion action for 2-flavor 
QCD with mps/my ~ 0.7. The result of the critical value of h at which the first order transition 
appears, h c , is about 0.06. The difference may be caused by the lattice discretization error due to 
small N t . We have defined the parameter h = 2Nf(lK\f) N < for the Wilson quark. Then, the critical 
point Khc corresponding h c decreases as K/ lc °< (h c /Nf) l ! Nt with Nf, and the truncation error from 
the higher order terms of the hopping parameter expansion in ffh becomes smaller as Nf increases. 

The remarkable point of this study is light quark mass dependence. In Fig. [|, we plot the 
results of the critical value h c as functions of (mps/my) 2 (left) and the quark mass defined by the 
PCAC relation m pcac (right). The slope of V e ff is computed fitting the data with 5 th or 6 th order 
polynomials. The open symbols are the results by 5 th order, and the filled symbols are those by 
6 th order. The difference is taken as the systematic error. The red lines are the upper and lower 
bounds of the critical h for the system with Nf = 16 massive quarks and no light quarks, which is 
determined by observing the hysteresis curves of the plaquette and Polyakov loop. It is found from 
these figures that the light quark mass dependence of the critical mass of heavy quarks is very small 
in the region we computed. The weak dependence suggests that the critical mass of heavy quark 
remains finite in the chiral limit of 2-flavors. Thus, the sign of the first order transition in 2-flavor 
QCD is not shown. 

4. Singularities in the complex chemical potential plane 

Let us turn to finite density QCD. The distribution function of 2-flavor QCD at finite density 
and the appearance of the double-well potential are discussed in Ref. [||], and the boundary of the 
first order transition region of (2 +Nf )-flavor QCD at finite density is computed in Ref. [|3]] using the 
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2-flavor QCD configurations generated with the p4-improved staggered quark action in Ref. 

The first order region is found to become larger as pi increases. 

In this section, we extend the potential analysis to the complex chemical potential, pi = pi^ e + 
ipii m . The numerical study of the singularities of 3f = 0 has a potential danger for large pi when 
use the reweighting method Jl2|]. We estimate the position of 3f = 0 from the probability 


we 


distribution of the complex phase, since 3? vanishes when the distribution function has two peaks or 


more and the contributions from these peaks cancel each other [|l 2\\. In this study, we investigate the 
plaquette distribution function and complex phase average with constraining the plaquette value. 
To avoid the sign problem, the Gaussian approximation [|8|] is applied 2 . 

We compute the probability distribution function of the plaquette P, 


w( 


>(P,P,H) = j@U 8{P-P) (det M) 2 e 6pNsiKp . 

We denote wo(P,fi) = w(P,j3,0). The normalized partition function is rewritten as 

np ,pi) 


(4.1) 


JT(J8,0) 


jfjjj J w(P,P,pi) dP = —py J R(P,pi)w 0 (P,P) dP. (4.2) 


Here, R(P,pi ) is the reweighting factor for finite pi defined by Eq. ( |2.3| ) with /3o = /3 and nij = °o for 
W-flavors. This R(P.pi) is independent of /J and R(P.pi) can be measured at any (5. In this study, 
all simulations are performed at pi= 0 and the effect of finite pi is introduced though the operator 
del M(pi)/ detM(O) measured on the configurations generated by the simulations at pi = 0. 

Because QCD has time-reflection symmetry, the partition function is invariant under a change 
from pi to —pi, i.e. R(P,—pi) = R(P.pi). Moreover, the quark determinant satisfies det M(—pi) = 
(det M(pi*))*. From these equations, we get 


[R(P,pi)]*=R(P,pi*). 


(4.3) 


This indicates that R(P,pi) is real valued function in the case of real pi, i.e. pi = pi*. Then, the 
plaquette probability distribution R(P,pi)wo(P, ft) is also real valued. However, once the imaginary 
part of pi becomes nonzero, R(P,pi ) is not a real number any more. We thus write the partition 
function, 

iTGM) = f e^ p ^\R(P,pi)\w 0 (Pl5) dP. (4.4) 

This complex phase (j) vanishes at pi\ m = 0, and 0 is monotonic function of pi\ m at small pi if we 
define the complex phase by a Taylor expansion of IndctA/(/r). 

If \R(P, pi)\wo(P,[5) becomes a double-peaked function which has two peaks of equal height at 
P+ and P , two phases coexist like first order phase transitions. When changing pi, the expectation 
value of plaquette changes discontinuously form P_e“? ip > to P\ e l0ip '^ at that point. And then, the 
partition function can be approximated by Z « (e'(’( p +) + e ‘<l > ( p -)) x (const.). When the difference 
between <j)(P+) and </>(P _) is equal to (2/i+ I )n with an integer n, the partition function will be 
vanishing, i.e. Lee-Yang zeros appear. The first-order-transition-like point in the complex plane 
corresponds to “Stokes line” in the infinite volume limit. We define the effective potential as 

^Preliminary results are presented in Ref. ||o|] 
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P 



P 


Figure 4: The derivative of the effective potential dV e s/dP at /j = 3.65 for (p^/T) 2 — (p lm /T) 2 = 1.0 
(left) and 4.0 (right). 


V e ff = \n \R(P. p ')\w()(P. /3)] and investigate the position of this would-be Stokes line in the complex 
p plane. 

We calculate these three quantities, |i?(P,ju)|, (f>(P.p) and wo(P,j3) by Monte-Carlo simula¬ 
tions. Because the exact calculation of the quark determinant is difficult except on small lattices, 
we estimate the quark determinant from the data of Taylor expansion coefficients of lndet M(p) 
up to 0(p 6 ) around p = 0 obtained by a simulation of 2-flavor QCD with p4-improved staggered 
quarks in Ref. [|l3]] . The truncation error has been discussed in Ref. [|J]. Notice that the complex 
phase 9 is not restricted to the range from — K to n because we define 6 by the Taylor expansion 
of IndetM. 

To compute Veff at finite p, we discuss the sign problem. We denote the quark determinant as 
N\\n\Ac{M(p)/ detM(O)] = F + i9. Histograms of 9 seem to be well-approximated by Gaussian 
functions. (See Fig. 1 (left) in Ref. [j|].) Here, we perform a cumulant expansion, 


(exp (F + i 9)) = (e F )exp 


i(9) - \({A9) 2 ) - j|-((A0) 3 } + i(AFA9) - l -(AF(A9) 2 ) + ■ 


(4.5) 


with AX =X- (X). 

For the case that the distribution of 9 is of Gaussian, the 0(9 n ) terms vanish for n > 2 in this 
equation. Moreover, since 9 ~ O(p) and F ~ 0(p 2 ), this expansion can be regarded as a power 
expansion in p. If the convergence of the expansion Eq. ([O]) is good, we can extract the phase 
factor e 1 ^ from R easily and the sign problem in \R\ is eliminated. We deal with the first two terms, 
i.e. i(9) and —((A0) 2 )/2, assuming the Gaussian distribution. The correlation terms between F 
and 9 are also neglected as a first step. Because we calculate the expectation value with fixed P and 
the values of F and P are strongly correlated, the A F may be small once P is fixed. Then, <j) ~ (9). 

We discuss dV^/dP instead of V e ff itself, again. dV^/dP at different /3 can be estimated by 
the equation, 

(P,Po) ~ 6(/3 — A))Afsite, (4.6) 

under the parameter change from /3 q to /3. 
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Figure 5: Left: The complex phase 0 for Ref/I/T) 2 = 1 as a function of P. Right: The blue lines are 
the would-be Stokes lines for each p. The values of p are shown near the blue lines. The red curve is the 
boundary where V e ff changes to double-well type in the complex fi/T plane. Below this line. Veil is always 
of single-well. 


If the effective potential Veff is a double-well function of P having two minima at P and P 
and one maximum at the middle Pq, dV e ff/dP is an S-shaped function and vanishes three times at 
P , P and Pq. The condition, 

[ P ° (dV efi /dP)dP = - [ P+ (dV e ff/ dP)dP, (4.7) 

7p_ Jp 0 

is satisfied when V e ff(P+) = V e ff(P_) at the transition point. In such a case, there exists a region of 
P where the derivative of dV^n/dP is negative. The results of the first derivative of V e ff(P,j8,ju) are 
shown in Fig. |] for various /I /T with Re[(fl/T) 2 } = (/i Re /T) 2 — (fi\ m /T) 2 = 1 (left) and 4 (right), 
where P = 3.65. We measured them at the peaks of the plaquette histograms and interpolated 
the data by a cubic spline method. In the region of large fli m /T, dV e s/dP becomes an S-shaped 
function. The result of 0 is plotted in Fig. [5] (left) for Rc\(ji/T) 2 ] = 1. It is a monotonically 
increasing function of P. 

We then investigate the boundary at which dV^j/dP changes to an S-shaped function from a 
monotonic function. The boundary is shown as a red line in Fig. |5] (right). Above this line, the 
region of P where d 2 V c ((/dP 2 = 0 appears. The constant part of dV^/dP is changed by varying /3. 
We then adjust /3 such that the depths of two minima of V e ff are the same using Eq. (p~6|). If the 
dashed lines in Fig. [| move to the horizontal axis by changing /3. Eq. (£0]) is satisfied. 

Contour plots of the /3 at which V e ff(/+) = V e rr(P .) are shown by blue curves in Fig. |5] (right). 
The values of p are indicated near the blue lines. Along the line for each /3, the double-well 
potential appears, and the contour curve is expected to turn into the Stokes line in the infinite 
volume limit. For finite volume, at points on that line, where the phase cancelation occurs, Lee- 
Yang zeros appear. This result indicates the existence of singularities in the region of large /./| m as 
well as the region of large ^ Re , and the boundary in the complex plane is closer to the origin /./ = 0 
than that (the square symbol) on the real axis. 
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5. Conclusions 

We studied the distribution function and effective potential of QCD with two light flavors and 
Nf massive flavors, aiming to understand phase structure of 2 and (2+l)-flavor QCD. Through the 
shape of the distribution function, we investigated the critical surface separating the first order 
transition and crossover regions in the parameter space of the light quark mass, heavy quark mass 
and the chemical potential. It is found that the critical mass becomes larger with Nf and the first 
order region becomes wider with the increasing chemical potential for large Nf. If (2+l)-flavor 
QCD has the same property, this gives the strong evidence for the existence of the critical point at 
finite density in the real world. 

The nature of the chiral phase transition in the 2-flavor massless limit is still open question. To 
study the chiral limit, we investigated the light quark mass dependence. If the transition is of first 
order, the critical vanishes before going to the 2-flavor massless limit. But, the critical jq, does 
not show such a behavior in the region we investigated. This implies that the critical mass of heavy 
quark remains finite even in the chiral limit of 2-flavors and 2-flavor QCD near the chiral limit is 
not in the first order transition region. 

We then discussed 2-flavor QCD with the complex /I by a numerical simulation. We found that 
there is a region where the plaquette distribution function has two peaks, suggesting the existence of 
singularities characterized by 2f = 0, in the region of large /q m as well as of large /Jr c . Combined 
with the analytical study of the complex chemical potential, the distribution of the singularities 
may provide important information about the QCD phase transition. 
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